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Introduction and background 


A digital image: An array of scalars or vectors. 

Scalar: Reflectance, temperature, range 
Vector: RGB, multispectral, hyperspectral 
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A digital image 



Landsat MSS image, 
courtesy of NASA 
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Image registration and 
image fusion 

Image registration is the process of spatially aligning 
two or more images of a scene. This spatial alignment 
is needed to fuse information in the images. 
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Landsat MSS Landsat TM 

Data courtesy of NASA 



Registered MSS & TM 
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Applications of image registration 

and image fusion 


Change detection 



Landsat 1 Landsat 2 Change image 

Data courtesy of NASA 
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Fusion of multimodal data 



Fused image 


Data courtesy of NASA 




Landsat TM bands 1 & 7 
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Image mosaicking 



Mosaicked image 


Two aerial images of Honolulu, HI. 
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Preprocessing Operations 
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Preprocessing operations 


All operations performed on images that 
improve the registration performance. These 
include: 

- Noise filtering 

- Deblurring 

- Region extraction 

- Edge detection 


A. Goshtasby 


Noise smoothing 


Given image f(x,y) and smoothing filter h(i,j), noise 
smoothing is defined by: 


The intensity of pixel (x,y) in the output is obtained 
from a weighted sum of intensities of pixels at and 
around (x,y) in the input. h(i,j) is the weight of pixel 
(x+ i,y+ j) in the neighborhood of (x,y), and the sum of 
the weights over all i and j is 1 . 


k 
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Mean filtering 

• When the weights defined by the filter are all the same, the 
operation is known as mean filtering. 

• Intensities of all pixels at and around a point in input have the same 
effect on the intensity at the same point in output. 

• This operation is not rotationally invariant if the filter kernel is not 
circular. 



Image containing Filter-radius 2 pixels Filter-radius 4 pixels 
Zero-mean noise Computed using FFT algorithm 
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Gaussian filtering 


If the weights in filter kernel represent 
Gaussian coefficients, a Gaussian filter is 
obtained. 

Gaussian filtering is effective when image noise is 
zero-mean. 

Gaussian filtering is rotationally invariant. 

A 2-D Gaussian can be decomposed into 2 1 -D 
Gaussians: G(x,y j = G(x) * G(y); therefore, filtering 
can be carried in 1 -D rather than in 2-D 



A. Goshtasby 
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Computation of mean and 
Gaussian filtering 

Although filtering is a convolution operation and can be computed using 
the FFT algorithm, since FFT considers an image is a periodic signal, if left 
and right image borders, or top and bottom image borders are not the same, 
artifacts will appear near the image borders. To avoid this, carry out the 
computations directly. 



Image containing Computed with FFT Computed directly 
Zero-mean noise Gaussian filter of o = 2 pixels 


A. Goshtasby 
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Image segmentation 


This is the process of partitioning an image 
into meaningful parts. There are two main 
approaches to image segmentation. 

- Methods that use information within regions 

- Methods that use information on the boundary 
between regions 


A. Goshtasby 


Intensity thresholding 


Assuming an image that contains objects with intensities that 
represent a Gaussian distribution and a background with 
intensities that represent a different Gaussian distribution, the 
objects can be separated from the background using the 
intensity at the valley between the two histogram modes. 

This method works well when an image contains 
homogeneous objects and a homogeneous background and the 
properties of the objects and the background are different. 



Original 


Smoothed 


Histogram 


Thresholded image 
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Intensity thresholding 



A Landsat image 



Thresholding at the 
valley between the 
first two peaks. 




average intensity of 
highest-gradient pixels. 
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Threshold selection 


The threshold value is computed by: 

- Finding the valley between the modes of the 
histogram of the image. 

- Finding the intensity that represents the average of 
intensities of high-gradient pixels. 

- Finding the intensity at which a change in the 
intensity will minimally change the segmentation 
result. 


A. Goshtasby 
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Edge detection 


Edge detection methods can be categorized into those 
that search for locally maximum image gradient 
magnitudes and those that search for zero-crossings 
of the Laplacian of an image. 

Methods that search for gradient peaks do not pick 
false edges but the ones picked could be 
disconnected. 

Methods that search of the zero-crossings of the 
Laplacian image find closed boundaries, but parts of 
the boundaries could be false. 


A. Goshtasby 
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LoG edge detector 


Determination of the LoG of an image 
involves computation of: 


LoG[f(x,y)' 


d 2 [f(x, y) * G (x,y)] d 2 [f(x,y ) * G(x, y)} 


= f{x,y)* 


dx 2 

d 2 G(x , y) 
<9x 2 


2 


+ f{x,y)* 


d 2 G(x,y) 


m2 


°r ? 


LoG[/(a:, ?/] ) = ^ ’ ★ G(i/) * f(x 9 y) + G(x) ★ 


Xv) 


* f ( x > 2/) 
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Edge detection by intensity ratio 






Zero-crossing edges can be 
determined by convolving 
the negative and positive 
parts of the LoG with an 
image separately and 
subtracting the convolved 
images and locating the zero- 
crossings. 

If instead of subtracting 
corresponding values in the 
convolved images, we divide 
the values and locate the 
one-crossings, intensity-ratio 
edges will be obtained. 



LoG 

operator in 
2-D 



Positive 

part 



Negative 

part 
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Edges in color images 


Edge detection in a color image can be considered edge 
detection in a 2-D vector field. 

If u(x,y) and v(x,y) represent color gradients in x and y 
directions, edges can be considered points where color 
gradients are locally maximum in the gradient direction. 

If R(x,y), G(x,y), and B(x,y) represent red, green, and blue 
color components at (x,y), respectively, color gradients are: 


u(x,y) 

v(x,y) 


dR(x, y) dG(x,y) 
r H ^ g 


dx 

9R(x, y) 

dy 


dx 

, 9G(x,y) 
r + 7 . g + 

dy 


dB(x,y) 

dx 

dB(x,y) 

dy 


r, g, and b are unit vectors along red, green, and blue axes, 
respectively, in the color space. 
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Color edges 


Gradient direction at (x,y) is the direction 
maximizing 


F(x,y ) = [u(a;, y) cos 0(x, y) + v(x, y) sin0(x, y)] 


2 


and is obtained from 


0(x, y) = 0.5 tan 


2u (x,y) ■ v(x, y) 


u(x, y) ■ u(x. y) - v(x, y) ■ v(x, y) 
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Edge detection in color images 



A color image 



Edges of the color image 
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Landmark Detection 

Arthur Goshtasby 
Wright State University and 
Image Registration and Fusion Systems 


Image features 


Points 


Lines 



Templates 
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Landmark detection 


They are 1 ) locally unique and 2) rotationally 
invariant. 
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Comers detected at different resolutions. 




Corners as landmarks 


Given image 1: 

1 . Compute image gradients in x and y directions: lx, Iy 

2. Compute square gradient matrix at each pixel (x,y): 



C( x,y) = 


IJ X 

Iyh 


771 


x y 

Vy 




where overbar implies average in a small neighborhood of 

(x,y). 

Compute eigenvalues of C(x,y) and if the smaller eigenvalue 
is locally maximum, take (x,y) as a comer. 


Examples 


Stable corners: Comers that persist over a wide 

range of scales/resolutions. 
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Spot centers as landmarks 


• Find points where 
Laplacian of Gaussian 
(LoG) of image is locally 
extremum. 

• Scale-invariant feature 
transform (SIFT): LoG 
points that are locally 
extremum in scale and 
space. 
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Locally unique and highly 
informative points as landmarks 

Given an image and a threshold distance d: 

1 . Find edges in the image. 

2. At each edge point compute entropy within a circular window 
centered at it. 

3. Remove edge points that have produced the lowest p% entropies. 

4. Compute the auto-correlation at the remaining edge points and 
save an edge point in list L, which is initially empty, if the largest 
auto-correlations there is locally minimum. 

5. Sort list L from the smallest to the largest auto-correlation. 

6. Consider the edge point on top of the list the next comer. 

7. Remove all edge points that are within distance d of the selected 
edge point from L. 

8. If there are still edge points in L, go to Step 6. Otherwise stop. 
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An Example 




Line intersections as landmarks 


Line detection 

• Least-squares 

• Hough transform 

• Image gradients 


Least-squares line fitting 


Given a set of points {( Xi,yt): i=l,...,n} find p and 0 of 
line p = xcos<9 +ysmd such that 

E = Y j (p-x i cos 0 - y t sin 0f 

i = 1 

is minimum. 



X 


X 


10 




An example 



Canny edges Detected lines 
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Hough transform 


A point (x,y) in the xy-space represents a unique 
line y = mx + b in the m^-space, and a point 
(m, b) in the mb-spacG represents to a unique line 
in the xy-space. 
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Hough transform 

Parameter m cannot be accurately determined for 
nearly vertical lines. So, instead of y=mx+b, the 
polar form of line p = xcos0 + ysinG is used. 


Algorithm: Given an edge image: 

1 . For each edge point (x,y), draw p = xcos0 + ysin0 
in a p0 array. As each sinusoidal is drawn, 
increment p0 entries in the array that fall on the 
sinusoid by 1 . 

2. Locate locally peak entries in the array. Each such 
entry in the array p0 detects a line in the image. 
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Line detection 
using image gradients 


1 . Determine the gradient magnitude and gradient 
direction of each pixel in the image. 

2. Group the pixels into regions with gradient 
directions in [-a, a] [a, 3a] . . , where a is a 
small angle, such as 5 degrees. 

3 . Remove regions that contain fewer than a 
required number of pixels. 

4. Fit a line to each remaining region by the 
weighted least-squares method, with the weights 
being the gradient magnitudes at the pixels. 
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An example 



Lines detected at two different values of a. 
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Region centers as landmarks 


•For some image types, regions are the best 
features to use in image registration. 

• Centers of gravity of regions are not sensitive to 
zero-mean noise. 

• Centers of gravity of regions can be determined 
with subpixel accuracy. 

• Centers of gravity of regions are affine invariant. 

• Regions can be obtained by various image 
segmentation methods. 
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An example 



Landsat MSS image 



Smoothed and thresholded image 
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Template centers as landmarks 


• Templates are image regions that are locally 
unique and informative. 

• They can be considered circular windows that 
are centered at detected landmarks. 

• Circular windows make the selected features 
independent of an image’s orientation. 
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An example 



Corners Templates 
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Similarity Measures 


Arthur Goshtasby 
Wright State University 
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Definitions 


A measure to quantify how well two patterns 
or two image windows match. 

Dissimilarity measure is also known as 
distance measure. 

A similarity/dissimilarity measure may or may 
not represent a metric. 

If SD is a metric similarity/dissimilarity, then 
SD(X,Y) = SD(Y,X), for windows X and Y. 


Pearson correlation coefficient 


Given intensities X = {xt : i = 1, ...,n} and 
Y = {yt : i = 1, 




r varies between + 1 and - 1 . + 1 shows perfect 
positive correlation and - 1 shows perfect 
negative correlation. 


Computation of correlation 


• Using Fourier transform 

C = T~ x [T{U)-T*{V)\ 



• Phase component of 


cross-power spectrum 



T{V)-F{V) 

: f(u)-P(v) 
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Spearman’s rho 


Suppose replacing intensities xt and yt with their ranks 
R(xi) and R(yi) and then calculating the Pearson 
correlation coefficient between the ranks. This is 
equivalent to finding: 

^ . 6Eti \R(xi)-R{yi) 

n(n 2 - 1) 

This measure is invariant to monotone 
intensity differences between images. 




Shannon’s mutual information 


Assuming pij is the probability that 
corresponding pixels in X and Y have 
intensities i and j, mutual information is 
defined by: 



z=0 j - 0 



Li and L2 norms 


Li norm is also known as sum of absolute 
intensity differences. 

n 

Li = I Xi - Vi 

i — 1 

L 2 norm is also known as the sum 
of squared intensity differences. 



Joint entropy 


• Having the joint probability density py of 
intensities i and j in X and Y, joint entropy is 
defined by: 


255 255 

D E = ~J2J2 Pij lo §2 Pij 

i=0 j = 0 

• De should be minimized to find correspondence. 
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Exclusive /-information 


Information exclusively contained in X and Y when 
observed jointly is known as /-information. 

Exclusive /information is related to joint entropy De 
and mutual information Smi by: 

D f (X,Y) =D E (X,Y) - S M i(X,Y) 



Properties 


Sensitivity to noise: Spearman’s p and Li norm. 

Sensitivity to intensity differences: Mutual 
information, joint entropy, and Spearman’s p. 

Sensitivity to blurring: Pearson’s correlation 
coefficient, Li norm, and L2 norm. 

Sensitivity to geometric difference: Pearson’s 
correlation coefficient and Li norm. 

Speed: Li norm, L2 norm, Pearson’s correlation 
coefficient. 


Related references 


• D. I. Bamea and H. F. Silverman, A class of algorithms for fast digital image 
registration, IEEE Trans. Computers, vol. 21, no. 2, pp. 179-186, 1972. 

• W. Kruskal, Ordinal measures of association, Journal of American Statistical 
Association, vol. 53, pp. 814-861, 1958. 

• C. D. Kuglin and D. C. Hines, The phase correlation image alignment method, Proc. 
Int'l Conf. Cybernetics and Society, pp. 163-165, 1975. 

• F. Maes, D. Vandermeulen, and P. Suetens, Medical image registration using mutual 
information, Proc. IEEE, vol. 91, no. 10, pp. 1699-1722, 2003. 

• K. Pearson, Contributions to the mathematical theory of evolution, III, Regression, 
heredity, and panmixia, Philosophical Transactions Royal Society London, Series A, 
vol. 187, pp. 253-318, 1896. 

• J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever, /^information measures in medical 
image registration, IEEE Trans. Medical Imaging, vol. 23, no. 12, pp. 1506-1518, 2004. 

• C. Spearman, The proof and measurement of association between two things, The 
American Journal of Psychology, vol. 15, no. l,pp. 72-101, 1904. 
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Image Features 


Arthur Goshtasby 
Wright State University 
Image Fusion Systems Research 


Types of image features 


Statistical 

Geometric 

Algebraic 

Frequency 

domain 


• Filter response 

• Differential 

• Fractal 
dimension 

• Information 
theoretic 


Statistical features 


Assuming p(i) is the probability that intensity 


appears in image: 
Mean: 

255 

/' = ip (®) 

i=0 

Skewness: 

OKK 

-j z o o 

7 = -3 S(* - /^) 3 P(*) 

i=0 


Variance: 

255 

V 2 = £(* - l4 2 p(i) 

i = 0 

Kurtosis: 

055 

« = ' ; 2(< - i<)'i>(i) - 3 
^ *=0 


Examples 




Image 


Feature obtained when using 
windows of radius 8 pixels. 




Skewness 


Kurtosis 
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Geometric features 


Invariant moments: 


Central moments: Assuming (xo,yo) is center of gravity of image f(x,y) 

m — 1 n — 1 

y Pq = Y Y^ x ~ x ° f{y - yo) q f(x , y) 

x=0 y = 0 

Invariant moments: (order 2) 


II 

Ii 


(H 20 + M02), 

(^20 - M 02 ) 2 + 4 / i ^ 


Examples 



Invariant moment Ii 


Invariant moment 


Algebraic features 


Given image f, if U is the column eigenvector 
system of f and V is the row eigenvector 
system of f, we can write: 


UfV 




£ = diag{(7\ , 02, . . . , a r ) , where 

1/2 

<7j = A/ is the zth largest singular value. 


Singular values 


• Rather than Oi, use ci/ci; it is invariant to image 
contrast. 


Example: 



( 72/(71 
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Frequency domain features 


Power spectrum features: If F(u,v) is the 
Fourier transform of image f(x,y), then 

( f>(u,v ) = F(u,v)F*(u,v) = \\F(u, t ;)|| 2 

are the power spectrum features of the image. 
Total power spectrum: 

J2J2\\ F ( u ’ v )\\ 2 



Filter responses 


• Consider 1-D filters: 


B 0 = 

Bi = 

B-2 = 

B-i = 

B a = 





and suppose creating 2-D filter Bij=Bi*Bj, then 
response of image to By can be used as a 
feature. 
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An example 


• Response to Bn: 
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Differential features 


Rotationally invariant differential operators 


Smoothed: /(x, y), 
Gradient Magnitude: 



2 

x 


+ /„(•'■• //)' 1 2 


Laplacian: f xx ( x , y ) + f yy ( X, y ) , 





i ' 

//kll 

1 - < 

9 1 





Smoothed 
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Spatial domain features 


Absolute center contrast: 


l 

M N - 1 


M- 1 N - 1 


\f( x iy) - fc 

;r=() y — 0 



Fractal dimension 


• It is a measure of roughness when considering 
intensities as height values. Change in surface 
area as a function of change in resolution is 
used as fractal dimension. 
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An example 



Fractal dimensions 
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Information theoretic features 


• Entropy: It measures information content in an 
image or window. 



E p (o lo s2 m 

i=0 







9 * 


♦ 


/ 
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Properties of features 


Invariance under blurring: Invariant moments, 
smoothed intensity 

Invariance under noise: Smoothed intensity, 
invariant moments, fractal dimension 

Invariance to change in contrast: Gradient 
filter response 

Invariance to rotation and scale: Smoothed 
intensity, invariant moments, Laplacian 


Related references 


• H. C. Andrews and C. L. Patterson, Singular value decomposition and digital image 
processing, IEEE Trans. Acoustics, Speech, and Signal Processing, pp. 26-53, 

1976 

• M. K. Hu, Visual pattern recognition by moment invariants, IEEE Trans. 
Information Theory, vol. 8, pp. 179-187, 1962. 

• J. J. Koenderink and A. J. van Doom, Representation of local geometry in the 
visual system, Biological Cybernetics, vol. 55, pp. 367-375, 1987. 

• K. I. Laws, Rapid texture identification, in Image Processing for Missile Guidance, 
Proc. SPIE , vol. 238, pp. 376-380, 1980. 

• A. Pentland, Fractal-based description of natural scenes, IEEE Trans. Pattern 
Analysis and Machine Intelligence, vol. 6, no. 6, pp. 661-674, 1984. 

• C. E. Shannon, The mathematical theory of Communication, in book with the same 
title by C. E. Shannon and W. Weaver, University of Illinois Press, Urbana, 1949, 
reprint 1998. 
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Image Descriptors 


Arthur Goshtasby 
Wright State University 
Image Fusion Systems Research 


What is an image descriptor? 


It is a feature vector. Traditionally, all 
features are taken to be of the same type, 
but this is a self-imposed restriction. 


• Histogram-based • Shape context 

• SIFT • Spin image 

• GLOH • RIFT 


Histogram-based descriptors 


A feature vector with components representing 
the histogram bin counts. 

Assuming Hi and H2 histograms of two 
windows, and the zth bin count in the 
histograms are Hi(i) and Hi(i), distance 
between the two histograms is determined 

from D(Hi,H2)—Ui mim{Hi(i),H2(i)}. 


Scale-invariant feature transform 

(SIFT) 

• Find SIFT points. 

• Estimate scale and orientation of neighborhood 
of each point. 

• Normalize a local neighborhood and subdivide 
into 4x4 blocks. 

• Create a gradient histogram for each block. 

• Concatenate the local histograms to create the 
descriptor. 


Finding SIFT points 


Convolve image with Gaussians of standard 
deviations V2, 2, 2 V2 , 4 ..., creating a stack of 
Gaussian smoothed images. 

Find difference of Gaussian (DoG) smoothed images 
in the stack to create another stack, approximating the 
LoG of the image at different resolutions. 

Consider the stack of the DoG images a volumetric 
images and find extremum points in the volume. 

Slice number of an extremum determines the scale of 
neighborhood of the point. 


Estimating orientation 


Knowing the slice number, n, of a point in the DoG 
stack, find the gradient magnitude and gradient 
direction of slice n in the Gaussian smoothed stack. 

Weigh gradient magnitudes by a Gaussian with 
standard deviation proportional to n. 

Find histogram of gradient directions within a square 
neighborhood of width proportional to n, using 
weighted gradient magnitudes as increments. 

Take the direction representing the histogram peak as 
the orientation of the point. 


Normalizing and subdividing 

neighborhood 



Gradient directions in each block are normalized with 
respect to the peak orientation and grouped into a 
Histogram with 8 bins. 




Gradient location and orientation 

histogram (GLOH) 


• Subdivide the neighborhood 
into a log-polar grid with 
radii 5, 11, 16, and 90- 
degree angular spacing. 

• Group gradient directions 
within each block into a 
histogram with 16 bins. 

• Overall, producing 272 
numbers, but reducing to 
128 by PC A. 
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Shape context 


Same as GLOH, but spacing 
the rings uniformly spaced in 
log space. This emphasizes 
the center part of a window. 

Five circular blocks are used 
and pixels between 
consecutive rings are 
grouped into 1 2 directions, 
producing overall 60 
numbers. 




Spin image 


Consider grouping intensities within each circular block 
into a histogram with 8 bins, and mapping the bin 
counts to column entries in a matrix. This will produce 
a rotation-invariant matrix known as spin image. 



Image 


Spin image 
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Rotation-invariant feature 
transform (RIFT) 

• If instead of image intensities, gradient 
directions that are measured radially are used, 
the obtained spin image has been named RIFT. 

• Note that gradient direction measured radially 
is independent of the orientation of the images, 
thus creating a rotation invariant descriptor. 
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How reliable are SIFT rotation 
and scale parameters? 


0 yBk T A 


Original Rotated by 30 degrees Scaled 


SIFT points with associating orietations and scales 


SIFT rotation parameters 

• Find MST of SIFT points before and after rotation. 

• Consider points with the same degree in the two MSTs. 

• Find difference in SIFT rotations at such points. 

• Find histogram of the rotational differences. 

• Peak in the histogram shows rotational difference between 
images. 

• SIFT rotation parameters 
appears useful in the 
absence of scale 
difference. 



SIFT rotation parameters 


• When the images have scale difference, 
rotation parameters are not very useful. 



0 30 359 


14 


SIFT scale parameters 


When the images are in the same scale, scale 
parameters are useful. Show below is the 
histogram of the scale ratio of MST points with 
the same degree in the original and rotated 
images. 



1.0 
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SIFT scale parameters 


• When the images have different scales, scale 
parameters are not useful. Shown below is the 
histogram of the scale ratios of original and 
scaled images. 



0.67 1.0 
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The Correspondence Problem 


Arthur Goshtasby 
Wright State University 
Image Registration and Fusion Systems 


Correspondence methods 


Point pattern matching: 

1) scene coherence 

2) clustering 

3) invariance 

Line matching 
Region matching: 

1) shape matching 

2) relaxation labeling 

3) chamfer matching 

Template matching: 

1) Similarity measures 

2) Coarse-to-fine methods 
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Point pattern matching 


Problem: Given two sets of points, 


determine the correspondence between them. 

V 

8 o ° 

- Information about only the locations of the 

Vo° 

points is available. 

°o %° 

- The points may contain positional error. 

- Some points may exist in only one of the sets. 


We will only consider the case where the two 

o 

o 

8° °° 

° o° 
°0 0 0 

point sets are related by the affine 

transformation. 

u o 




Point pattern matching using 

scene coherence 


If three points in the two sets are aligned, because of the 
scene coherence, the remaining points in the two sets will 
also align. 

RANSAC Algorithm: Find three corresponding points. 

1. Take 3 points from set 1 and 3 points from set 2 
and find the affine transformation that aligns them. 

2. Count the number of other points in the two sets 
that also align with the obtained transformation. 

3. If the count is sufficiently high, stop the process. 
Otherwise, go to step 1. 




To speed up the search, limit the point combinations to 
those falling on the convex-hulls or the minimum-spanning 
trees of the two sets. 
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Point pattern matching using 

clustering 


Algorithm: 

1 . Create accumulators a[ ] — f[ ] and initialize the 
entries to 0. 

2. From point triples in the two sets calculate 
parameters a-f of 

X=ax+by+c, 

Y=dx+ey+f. 

3. Increment entries a-f of accumulators a[ ] - f[ ], 
respectively, by 1 . 

4. Repeat Steps 2 and 3 a sufficiently large number of 
times. 

5. Locate the entry with the highest counts in a[ ] ... 
f[ ]. They show parameters a-f of the registration. 




f 
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Point pattern matching using 

affine invariance 


If two point sets are related by an affine 
transformation, knowing three corresponding 
points in the two sets (pi,qi), (p2,q2), (p3,q3), 
the relation between corresponding points 
( p,q ) in the sets can be written as 

p = pi + ai(pi - pi) + a2(p3 - pi) 
q= qi+ ai(q 2 - qi) + a2(q3 - qi) 






Point pattern matching using 
affine invariance 


Algorithm: 

1 . Create two 2-D accumulator arrays Hi [ ] 
and H 2 [ ]. 

2. Select three points in set 1 and for each 
additional point in set 1 calculate ai and a .2 
and increment entry [a,,a 2 \ of Hi [] by 1. 

3 . Select three points in set 2 and for each 
additional point in set 2 calculate ai and a 2 
and increment entry [a,,a 2 \ of H 2 [ ] by 1. 

4. Find the similarity between Hi and H 2 . If 
the similarity is sufficiently high, take the 
point triples selected in Steps 2 and 3 as 
corresponding points and stop. Otherwise, 
go to Step 2. 
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Line matching 


Assumption: Images are related by a 
rigid transformation (unknown 
translation and rotation). 

Algorithm: 

1 . Find the rotational difference 
between the line sets. 

2. Correct the orientation of set 2 with 
respect to set 1 . 

3. Find the translational difference 
between the line sets. 
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Region matching 


• Shape matching 

- Fourier descriptors 

- Invariant moments 

- Shape matrices 

- Relaxation labeling 

- Chamfer matching 

• Distance transform 




9 


Fourier descriptors 


Given pixels on the boundary of a region 
{(xi,yi): i=0, the xy coordinates of the 

pixels along the boundary can be considered 
samples from a periodic signal. Letting 
Zi-x&jyi, where j= V^l, Fourier descriptor is 
computed from 

1 N - 1 

G =— Z^ ex Vi-finkUN), k=0,...,N-l 

N i=0 


Invariant moments 


Given boundary pixels {(xuyi): i=0,...,N-l }, the (p+q) th 
order moment is defined by 


N - 1 


m 


pq 


= 2>,W, 


i = 0 


The (p+ q ) th order central moment of the boundary is 
defined by N _ } 

:=^(Xi-x) p (yi-y) q fi 


u 


Invariant moments: 

Cli = U20 U02 


i=0 


0,7 — (3ll2i — U 03 ) ... 


Shape matrix 


• Resampling a shape into a representation that 
is independent of its position, orientation, and 
scale. 


0 12 3 4 


/// >r"" 77\ \ \ 0 

11111 

((1 f \~J>— X. -L \ 1 

110 10 

1 | / f \qt\i 2] 3) to 2 

11110 

\ U \ vv 7 / / 3 

11110 

\ \\ L/i 4 

11110 


1110 0 
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Relaxation labeling 


• Denote regions in the reference image by 
{cn:i=l f and regions in the sensed image by 
{bi:i=l,...,n}. 

• If af s denote objects and bi s denote labels, the 
correspondence problem becomes that of finding 
labels for the objects such that a compatibility 
condition is satisfied. 

• Initially assign labels to the objects with probabilities 
proportional to their similarities. 

Pi(bj): similarity between regions at and bj. 

• Iteratively revise the label probabilities until they 
converge to either 0 or 1 . 


C 

T 
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13 

4 
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Chamfer matching 


Given two binary images with translational 
differences : 

1 . Initially position the sensed image within 
the reference image at (i,j) based on some 
available information. 

2. Determine the similarity between the two 
images using distances of closest object 
points. 

3. Reposition the sensed image within the 
reference image at the eight neighbors of 
(i,j) and determine the image similarities at 
these eight positions. 

4. If highest similarity is obtained when sensed 
image is at (i,j), stop. Otherwise, move the 
sensed image to the neighbor of (i,j) that 
produces the highest similarity and go to 
Step 2. 




Sensed 


A 


Initialization of 
sensed in reference 
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Distance transform 


• Given a binary image, assign to a 
background pixel a value proportional 
to its distance to the object point closest 
to it. 

• To speed up the computations, integer 
distances may be used, but note that 
integer distances involve errors and 
make distances dependent on the 
orientation of the image. 



Distance transform of a 
single point 



Isovalued distances 
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Distance transform 


To make distances 
rotationally invariant, use 
actual Euclidean distances. 

Distance transform in its 
current form is very 
sensitive to noise. 




Euclidean distance transform Eulidean distance transform 


of a circle of a circle and 5 noisy points 
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Instead of saving a 
single distance at a pixel, 
save a weighted sum of 
distances. 

This can be implemented 
by smoothing the binary 
image with a Gaussian 
and inverting the values. 



Circle Circle + 5 points 



points 5 points 
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Template matching 


This is same as chamfer matching except that 
the images are no longer binary. 

Similarity measures 

- Sum of absolute differences 

- Cross-correlation coefficient 

- Mutual information 

Gaussian weighted templates 
Coarse-to-fine approaches 


Similarity measures 


• Given two gray scale images with translational 
differences, shift one image over the other and at each 
shift position determine the similarity between the two. 


Sum of absolute differences 

m — 1 n — 1 

D(x, y) = ^2 I ~ fw[i + xj + y] 

i = 0 j = 0 

Cross-correlation coefficient (p = /- /mean) 

EEt/ E"=o 9t [i, j]ffw [i + X, j + y\ 


S{x,y) = 


{EEo 1 e;=o & 2 m EEo 1 e u + x -j + »]} 


1/2 
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Mutual information 


Create a 2-D histogram with entry {a, b\ 
showing the number of pixels in the template 
with intensity a that align with intensity b in the 
window. 

Divide the counts by the number of pixels in 
template to obtain joint probabilities. 

Then, compute: 


y 


Template 


Window 


y 

u 

\Cl 

\U 


2-D 

histogram 


255 255 


/(£, ui) 



Ptw{a,b ) log 


a =0 6=0 


Pf w &) 

P,(a)PJb) 
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Gaussian weighted templates 


• Instead of treating all pixels in a 
template similarly, give higher weights 
to pixels that are closer to the template 
center. 

• This makes the process less dependent 
on image orientation when using rectangular templates. 

• It also reduces the effect of geometric difference between 
images. 
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On template shape and size 


• Always take circular templates. 
This will make the process less 
dependent on image orientation. 


• Set template size proportional to the 
information content in the template. 
A smaller template size is sufficient 
when the template is highly detailed 
compared to when it covers a rather 
homogeneous area. 



s/ x 



For a detailed area 



For a less detailed area 
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Coarse-to-fine matching 


• At image level: Reduce the size of images, 
find the correspondences, and determine 
the transformation parameters. Use the 
transformation to resample the images at 
one level higher resolution. Repeat the 
process until images at the highest 
resolution are registered. 

• At template level: Either use smaller 
templates or cheaper similarity measures to 
find candidate match positions. Then find 
the best match position from among the 
candidates using a larger template and/or a 
more expensive similarity measure. 


< I 

i f ! 1 



Image pyramid 
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• During search: Use large steps to find 
possible match positions. Then use finer 
steps in the neighborhood of the likely 
matches to refine the final match position. 

• Use partial information: If a large set of 
landmarks is given, first use a subset of 
the landmarks to find approximate 
registration parameters and then verify the 
correctness of the registration and refine 
the parameters using all landmarks. 


Feature matching references 

1. A. Goshtasby and G. C. Stockman, Point pattern matching using convex-hull edges, IEEE 
Transactions on Systems, Man, and Cybernetics, 15(5):63 1-637 (1985). 

2. W. J. Rucklidge, Efficiently locating objects using the Hausdorff distance, Int'l J. 
Computer Vision , 24(3):25 1-270 (1997). 

3. G. Stockman, S. Kopstein, and S. Benett, Matching images to models for registration and 
object detection via clustering, IEEE Trans. Pattern Analysis and Machine Intelligence, 
4(3):229-241 (1982). 

4. J. L. Mundy and A. Zisserman, Geometric Invariance in Computer Vision, The MIT Press, 
Cambridge, MA (1992). 

5. J. Kittler and J. Illingworth, Relaxation labeling algorithms - A review, Image and Vision 
Computing, 3 (4):206-216 (1985). 

6. L. S. Shapiro and J. M. Brady, Feature-based correspondence: An eigenvector approach, 
Image and Vision Computing, 10(5): 283-288, 1992. 

7. A. Goshtasby, S. H. Gage, and J. F. Bartholic, A two-stage cross correlation approach to 
template matching, IEEE Transactions on Pattern Analysis and Machine Intelligence, 

6(3) :3 74-3 78 (1984). 

8. D. Holtkamp and A. Goshtasby, Precision registration and mosaicking of multicamera 
images, IEEE Trans. Geoscience and Remote Sensing, 47(10):3446-3455, 2009. 
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Transformation Functions 


Arthur Goshtasby 
Wright State University 
Image Registration and Fusion Systems 


Problem description 


Given the coordinates of a set of 
corresponding points in the images: 

{ (xi,yt) (Xi,Yi ): i=l,...,N) 

we want to determine function f(x,y) 
with components fx(x,y) and fy(x,y) 
such that 

Xi = fx(Xi,yi), 

Yi = fy(xi,yi) , i = 




(X,V 




Approach 


Rearrange the coordinates of corresponding 
points into two sets of 3-D points: 

{(xi,yi,Xi): i=l,...,N }, 

{(xi,yi,Yi): i=l,...,N}, 

then, fx and fy can be considered two single- 
valued surfaces fitting to two sets of 3-D 
points. 

We will consider the problem of finding 
function f(x,y) that approximates/interpolates 

{(xi,yi,fi): i=l,...,N }. 





Transformation functions 


• Translation 

• Rigid 

• Similarity 

• Affine 

• Projective 


• Thin-plate spline 

• Multiquadric 

• Weighted mean 

• Piecewise methods 

• Weighted linear 


Translation 


• Corresponding image 
points are related by 


X= x + h 
Y=y + k 


• One pair of corresponding 
points is sufficient to 
determine the registration 
parameters. 
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Rigid transformation 






X= x cos(0) —y sin(0) + h 

Y= x sin(0) +y cos(0) + k 

0 and (h,k) are the rotational and translational 
differences between the images. 

Knowing minimum of two corresponding points 
in the images these parameters can be 
determined. 

This transformation is useful when registering 
images as rigid bodies. 

Under rigid transformation shape and size are 
unchanged. 



Reference 



Sensed 
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Similarity transformation 


X= s x cos(0) — sy sin(0) + h 
Y= s x sin(0) + 57 cos(0) + k 

• s, 0 , and (h, k) are the scaling, rotational, and 
translational differences between the images. 

• Minimum two corresponding points in the 
images are required to determine s, 6, h, and k. 

• Angles are preserved under the similarity 
transformation. 

• This transformation is useful when registering 
distant orthographic images of flat scenes. 



Reference 



Sensed 
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Rigid registration example 



Registration of MR and PET images of the same person. 
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Affine transformation 


X= ax+by+c 
Y= dx+ ey + f 

• This transformation is a combination of shearing and 
similarity transformations. 

• The components of an affine transformation depend on 
each other. 

• Under the affine transformation parallel lines remain 
parallel. 

• When the components are made independent, the 
transformation becomes a linear transformation. 

• Knowing the coordinates of three non-colinear 
corresponding points in the images the six parameters of 
the transformation can be determined. 

• This transformation is useful when registering images 
taken from a distant platform of a flat scene. 



Reference 




Sensed 
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Affine registration example 



Registered 

lunar 

images. 




Subtracted 

registered 

lunar 

images. 
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Projective transformation 


X— (ax+by+c)/(dx+ey+ 1) 

Y= (fx+ gy+h)/(dx+ ey+ 1) 

• Knowing the coordinates of four non-colinear 
corresponding points in the images, parameters 
a- h can be determined. 

• This transformation is useful when registering 
images obtained from different views of a flat 
scene. 

• Under the projective transformation straight lines 
remain straight. 

• If images are from camera very far from a flat 
scene, projective transformation can be replaced by 
the affine transformation when registering the 
images. 




Reference 




Sensed 
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Thin-plate spline (TPS) 


N 

f(x, y) = A 1 + A 2 x + A 3 y + ^ Ftf In r 2 

i — 1 

where = (x-x ; ) 2 + (.y-.y,) 2 + c / 2 
• Parameters ^ 7 , ^ 2 , A3, and Ft: i=l, ...N are 

determined using {(xi,yi,fi): i=l,...,N} and the 
following three constraints: 


i:^ = 0 

r,^=o 



Sensed 
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Properties of TPS 


• Logarithmic functions are radially symmetric, so if 
the control point are not uniformly spaced, large 
errors may be obtained away from the control points. 

• TPS is useful when 

- the local geometric difference between images is not large 

- the control points are rather uniformly spaced 

- the density of the control points does not change 

- the number of corresponding control points is not very 
high. 
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Examples 




Sensed 

image. 



Registered 
using 
random 
grid points. 
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Multiquadric 


N 


Radial basis functions: 


When multiquadrics: 

When inverse multiquadrics: 

Note that whenGaussian: 


Ax,y)=YFA(*>y) 


i - 1 


Ri(*>y) - f( x ~ x i) 2 + (y~yt) 2 + d 2 ] lf 2 

Ri( x >y) = [( x ~ x i) 2 + (y~yi) 2 + d 2 ] 1/2 

Rt(x.y) = exp{- (x ~ Xi)2+( y~ yi)2 } 

2(7 


Multiquadric basis functions are also radially symmetric, therefore, their 
properties are similar to those of TPS. 
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Examples 


MAX: 0.0; RMS: 0.0 



Registration using MQ basis functions. 
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Weighted mean methods 


f(x, y) = X fPi (*> y) 

i = 1 

where b t ( x , y) = R t (x, y)/Yj, I, ^ (*> >’) 

and Ri(x,y) is a monotonically decreasing radial 
function. 

This is an approximation method; therefore, there 
no need to solve a system of equations. 


Problem with weighted mean 

When the weight functions are narrow, the functions produce 
flat spots at and in the neighborhood of the data points. 

Flat spots imply that many points in the sensed image map to 
the same point in the reference image, resulting in large 
registration errors. 

Consider points: (0,0,0); (0,1,0); (1,1,0); (1,0,0); (0.5, 0.5, 0.5) 



Narrow width 


Wider width 


Density of points 
(narrow width) 


Density of points 
(wider width) 
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Remedy 


• To produce a uniform density, use parametric 
surfaces as components of the 
transformation. 


X=fi(u,v) 

( 1 ) 

X= f2(u,v) 

( 2 ) 

y = fs(u,v) 

( 3 ) 


• For a given (x ,y), find corresponding (u,v) 
from (2) and (3). Then, use (u,v) in (1) to 
findX. 



1 



u 
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Resampled 

using 

single-valued 

surfaces. 


Resampled 

using 

parametric 

surfaces. 
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Examples 


MAX: 1.4; RMS: 0.6 MAX: 9.0; RMS: 1.3 




Registered using the weighted mean method. 
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Piecewise linear and 
piecewise cubic functions 

Algorithm: 

1 . Triangulate the points in the reference 
image. 

2. From the point correspondences, find 
corresponding triangles in the sensed 
image. 

3 . Map triangular regions in the sensed 
image one by one to the corresponding 
triangular regions in the sensed image 
using linear or cubic functions. 
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Examples 


MAX: 0.0; RMS: 0.0 
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Registered using piecewise linear functions. 
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Weighted linear method 


Note that weighted mean is a weighted sum of constants. 

Instead of using a weighted sum of constants, use a weighted 
sum of linear functions, each representing the geometric 
difference between the images locally. 

f(*> y)=Yjfi ( * ’ y) b t ( x > y) 

1=1 

wher e f. (x, y) = a t x + b t y + c t 

The local functions at (xi,yi) is determined by fitting a plane 
to (x i,yi,fi) and at least two other points in vicinity of (xi,yi) . 


Properties of weighted linear 

functions 


• The weighted mean method produces horizontal 
spots because a surface is obtained from a 
weighted sum of constants (horizontal planes). 

• In weighted linear, since the planes can have any 
orientation, horizontal spots are not obtained. 
This implies parametric surfaces are not needed 
to find the components of a transformation. 

• If rational weights are used, because the weights 
stretch toward the gaps, the functions can adjust 
themselves to the irregular spacing of the points. 

• If the width of the weight functions are set 
proportional to the local density of points, the 
functions can adjust themselves to the local 
density of the points also. 



Weighted mean 
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Examples 



MAX: 3.6; RMS: 0.7 
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Computational complexity 


Type of Transformation 

Computational Complexity 

Similarity 

0(N) + 0(n 2 ) 

Affine and Projective 

0(N) + 0(n 2 ) 

Thin-Plate Splines (TPS) 

0(N 2 ) + 0(n 2 ) 

Multiquadrics (MQ) 

0(N 2 ) + 0(n 2 ) 

Piecewise Methods 

0(NlogN) + 0(n 2 ) 

Weighted Mean 

0(n 2 N) 

Weighted Linear 

0(n 2 N) 
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What do transformation functions 
tell us about the images? 


• They contain information about 
geometric difference between 
images. 

fx(x,y) = x - 8sin(y/l 6) 
fy(x,y) = y +4cos(x/32 ) 


• They predict the geometry of 
the scene. 


Plot of fx(x,y) - x. 



Plot of f y (x,y) -y. 
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Detecting the mismatches 


Transformation functions contain information about the 
mismatches. 


> 

> 

) 

I 

Plot of fx(x,y j - x, using 5 
incorrect correspondences 


i 







Plot of fy(x,y) -y, using 5 
incorrect correspondences 


Generating image flow 

• Transformation functions can be used to generate flow diagrams, 
depicting the local motion or deformation of the scene 
from one image to the next. 




Volumetric registration 


Image flow 
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Image Resampling 


Arthur Goshtasby 
Wright State University 
Image Fusion Systems Research 


Problem description 

Problem: Knowing image intensities at discrete coordinates, 
determine intensities at arbitrary coordinates. 

Approach: Fit a surface to intensities at discrete coordinates and 
estimate the surface value at desired coordinates. 




Resampling methods 


• Nearest neighbor 

• Bilinear interpolation 

• Cubic convolution 

• Cubic spline interpolation 

• Radially symmetric kernels 


Nearest-neighbor resampling 


Set the intensity at (X, Y) to the intensity of the pixel 
closest to it: [round(X), round(F)]. 


Example: 



This method is very fast, but it produces aliasing 
effects along edges. 


Original 


Resampled 


Bilinear interpolation 

Assuming u and v are integer parts of X and Y, 
respectively, bilinear interpolation is defined by 

I(X, Y)= Wu,vl(u,v)+ Wu+i,vl(u+ 1, v) 

+ Wt/.vi ll(u,v+ 1 )+ Wu+l,v±ll(u+ l,v+ 1) 

where 

Wu,^(u+1-X)(v+1-Y) 

Wu+i,v=(X-u)(v+ 1-Y) 

Wu,v+i= (u+ l-X)(Y-v) 

Wu+i,v+i=(X-u)(Y-v) 




Bilinear example 



Original Resampled 
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Cubic convolution 


Cubic convolution can be computed row-by-row and 
then column-by-column. 

Assuming intensities at u-1, u, u+1, u+2 are I(u-l), 
I(u), I(u+1), I(u+2), intensity at A is estimated from 
f(X)=I(u-l)f,+I(u)f+I(u+ l)f+I(u+2)f 


where 


/-i 

/ o 
fi 
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1 +t 1 

2 2 

3 o 5 o , 

—t s tr -I- 1 

2 2 
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» 1 

r ( U+ 1 pV) 

1 It 


PtlMJ 


(u+Zv) 



Cubic convolution example 



Original 


Resampled 
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Comparison 


Suppose an image is rotated with a 10-degree increment 36 
times and then subtracted from the original image. 

Only nearest-neighbor preserves original image intensities. 
Bilinear and cubic convolution change image intensities. 



Nearest-neighbor 


Bilinear 


Cubic convolution 



Cubic spline interpolation 


• In cubic spline also, computation is carried out 
row-by-row and then column-by-column. 

• Assuming intensities at i = -1, 0, 1, 2 are 1 1, Io, 
h, L, the cubic B-spline curve approximating 
the intensities in the range 0 < u < 1 is 

/(«)=!>,(») 

i=-l 
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where b.^u) = (-u 3 + 3u 2 - 3u + l)/6, 
bo(u) = (3 u 3 — 6u 2 + 4) / 6, 

6i(u) = (— 3u 3 + 3 m 2 + 3u + l)/6, 

b-2{u) — u 3 /6. 

For the curve to interpolate the intensities, it is 
required to determine new intensities I'j s using Ls. 

j =- 1 

Instead of solving a very large system of equations, 
inverse filtering may be used to find the I'js. 


Radially symmetric kernels 

• Cubic convolution and cubic spline methods 
are not rotationally invariant. To obtain a 
rotationally invariant resampling, radially 
symmetric kernels are needed. 

• An example of a radially symmetric kernel: 


fin) = { 


1 — Sr‘f + 2 rf, 0 < r t < 1, 

0 . ri> 1 , 
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The influence of a pixel on interpolating function is 
radially symmetric. 

The resampled value at a point is obtained from the 
values of pixels within a circular area centered at the 
point. 




Computational complexity 


Type of Resampling 

Computational Complexity 

Nearest-Neighbor 

0(n 2 ) 

Bilinear Interpolation 

0{n 2 ) 

Cubic Convolution 

0(n 2 ) 

Cubic Spline, Direct Computation 

0(n 4 ) 

Cubic Spline, Using FFT 

0(n 3 log n) 

Radial Functions with Local Support 

0(n 4 ) 

Gaussian, Using FFT 

0(n 2 log n) 


For an n*n image. 
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Need for Fast and Accurate 
Image Registration 

Earth Science studies, e.g.: 

- Predicting crop yield 

- Evaluating climate change over multiple scales 

- Locating arable land and water resources 

- Monitoring pollution 

- Understanding the impact of human activity on major 
Earth ecosystems, etc. 

Global and repetitive measurements from a wide 
variety of satellite remote sensing systems 


Some Examples of Complementary Earth Science Missions 


Instrument 
(Spat. ResoL) 

Number of 
Channels 

AVHRR (D) 

(1.1 km) 

5 Channels 

TRMM/VIRS 

(2 km) 

5 Channels 

Landsat4-MSS 4 Channels 

(80 m) 

Landsat5&7-TM&ETM+ 

(30 m) 

7 Channels 

L ands at7-Panchromatic 

(15m) 

IRS-1 

LISS-I (73m) - 

4 Channels 
LISS-2 (36.5m) 

JERS-1 

8 Channels 


(Chl-4:18m; Ch5-8:24m) 


SPOT-HRV Panchromatic 
(10m) 1 Channel 

Spot-HRV Multispectral 

(20 m) 3 Channels 


MODIS 36 Channels 

(Chl-2:250 m;3-7:500m;8-36:lkm) 


EO/1 

ALI-MultiSpectr. 9 Channels (30m) 

ALI-Panchrom. 1 Channel (10m) 

Hyperion 220 Channels 

(30m) 

LAC 256 Channels 

(250m) 


IKONOS-Panchromatic 

(lm) 1 Channel 

IKONOS-MS 4 Channels (4m) 
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i i 


ASTER 14 Channels 

(Chl-3:15m;4-9:30m;10-14:90m) 


czcs 

Ofan) 


SeaWiFS (D) 
(1.1 km) 


8 Channels 


I I 


3 


4 


5-9 


10,11 

12 


13,14 


i i ~r 


i i 


i i 


i i 


TOVS-HmS2 (D) 20 Channels 

(15 km) 


I I 


I i 


GOES 5 Channels 

(1 km:l, 4km:2,4&5. 8km:3) 


METEOSAT 3 Channels 

(V:2.5km,WV&IR:5km) 


Visible 


Water 

Vapor 


IR 



Landsat ETM and IKONOS Registration 

US, Virginia Coast 



IKONOS 
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Image Processing Framework for Remotely Sensed Data 


Absorption 
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Feature 
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Signal to Noise at Wavelength X : 
(S/N^D^H/YT^L* 

Where 

D^: detectivity (measures detector performance quality) 

|i instantaneous field of view 
11 flying height of the spacecraft 
V : velocity of Ihe spacecraft 

\ : spectral bandwidth of the channel (spectral resolution) 
1^ : spectral radiance of pound feature 

= > TfibimU ftrm tYfi mml ami smual raitlutim, ox- 

To maintain the same SNR while improving spatial resolution fri' 
a factor of 4 ft.c., tier minus fi b a factor of If, we mutt degrade 
the spectral resolution hy a factor off (i.e.. increase Aft- a factor of 4) 
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The role of Image Registration 
in the Processing of Remotely Sensed Data 


Essential for spatial and radiometric calibration of 
multitemporal measurements for creating long-term 
phenomenon tracking data 

Used for accurate change detection: 

(Towsnhend et al, 1992) and (Dai & Khorram, 1998): small error in 
registration may have a large impact on global change measurements 
accuracy 

* 

- e.g., 1 pixel misregistration error => 50% error in NDVI 

computation (using 250m MODIS data) 

Basis for extrapolating data throughout several 
scales for multi-scale phenomena (distinguish 
between natural and human-induced) 


Normalized Difference Vegetation Index 


Classifying Image Registration Utilization 

Multimodal registration , for integrating 
complementary information from multiple sensors 


Multitemporal registration, for change detection and 
Earth resource surveying 

Viewpoint registration, for landmark navigation, 
formation flying (sensor web) and planet exploration 


Template registration, for content-based searching 
or map updating 
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Image Registration Requirements 


High Accuracy: Goal of sub-pixel accuracy 
Consistency: Robustness to recurring use 

Speed and High-Level of Autonomy: Needed for 

- Large amounts of data 

- Near- or Near-real time applications (e.g., disaster 
management) 


Systematic and Precision Corrections 


Navigation or Model-Based Systematic Correction 

- Orbital, attitude, platform/sensor geometric relationship, 
sensor characteristics, Earth model, etc. 

Image Registration or Feature-Based Precision 
Correction 

- Navigation within a few pixels accuracy 

- Image registration using selected features (or Control 
Points) to refine geo-location accuracy 

Two approaches 

1 . Image registration as post-processing 

2. Navigation and image registration in a closed loop 


Systematic and Precision Corrections 

AVHRR Example 
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Systematic and Precision Corrections 
AVHRR Example (cont.) 







Challenges in Registration 
of Remotely Sensed Data 

Image registration developed in other 
domains (medical, military, etc.) not always 
applicable 

- Variety in the types of sensor data and the 
conditions of data acquisition 

- Size of the data 

- Lack of a known image model 

- Lack of well-distributed “fiducial point” resulting 
in the difficulty to validate image registration 
methods in the remote sensing domain 

» use synthetic data, “ground truth", finer resolution dat|a 
and “circular" registrations 


Other Challenges Facing Image Registration 
In the Remote Sensing Domain 

• Navigation error 

- Historical satellites (e.g., Landsat-5 compared to 
Landsat-7) 

- Following a maneuver (e.g., star tracking) 

- Need for sub-pixel accuracy 

• Atmospheric and cloud interactions 

• Multitemporal effects 

• Terrain/relief effect 

• Multisensor data with different spatial and 
spectral resolutions 


Atmospheric and Cloud Interactions 

Baja Peninsula, California; 4 different times of the day (GOES-8) 

(Reproduced from Le Moigne & Eastman, 2005) 
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Multitemporal Effects 

Mississippi and Ohio Rivers before & after Flood of Spring 2002 (Terra/MODIS) 

(Reproduced from Le Moigne & Eastman, 2005) 



April 25, 2D02 


May 13.2C02 
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Relief Effect 

SAR and Landsat-TM Data of Lope Area, Gabon, Africa 

(Reproduced from Le Moigne et al., 2001) 
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Precision Correction in Operational Systems 


• Operational Environment 

• Platform/sensor models integrated 

• Historical data available for statistics/modeling 

• Robustness and consistency over time is a requirement 

• Operational Needs 

• Systematic correction (close to 1 pixel) using navigation model 

• Precision correction (less than 1 pixel) used to: 

- Check navigation model and ephemeris data 

- Perform band to band geometric calibration 

- Perform radiometric calibration of new sensor (relative to old one) 

• General Characteristics 

• Use database of Ground Control Points (GCP) or Chips 

• Normalized Cross-Correlation (NCC) is the most common similarity measure 

• Digital Elevation Model (DEM) is rarely integrated in the registration process 

• Cloud masking usually integrated 

• Errors in the [0.15-0.5] range 
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Precision Correction in Operational Systems 

Some Examples - Highlights 

• AVHRR: AUTONAV algorithm computes attitude corrections using Maximum Cross-Correlation 
(MCC) method between sequential images 

• GOES/METEOSAT: CPs and NOAA Shoreline database (GSHHS) used to match edges 
extracted from meteorological images 

• LANDSAT: CP image chips (1m orthorectified) using Gaussian pyramid, automatic Moravec 
window extraction and NCC or Mutual Information 

• MISR: Database of 120 GCPs (each a collection of nine geolocated image patches of a well- 
defined and easily identifiable ground features, from Landsat, terrain-corrected, data) &ray casting 
simulation software 

• MODIS: Biases and trends in the sensor orientation determined from automated control point (CP) 
matching and removed by updating models of the spacecraft and instrument orientation; finer 
CGPs from Landsat TM and ETM aggregated using PSFs and correlated with NCC 

• SEAWIFS: Reference catalog of islands GCPs and matching using spectral classification and 
clustering of data, “nearest neighbor” and pattern matching techniques 

• SPOT: Reference3D™ using DEM ortho-rectified simulated reference image in focal plane 
geometry, matching of input image to simulated using NCC and resampling into a cartographic 
reference frame 

• VEGETATION: Database of CPs from SPOT for VEGETATION1 and VEGETATION1 for 
VEGETATION2; Matching by NCC 
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Image Registration at NASA GSFC 


Operational Environment 

• Platform/sensor models integrated 

• Historical data available for statistics/modeling 

• Robustness and consistency over time is a requirement 

Operational Needs 

• Systematic correction (close to 1 pixel) using navigation model 

• Precision correction (less than 1 pixel) used to: 

- Check navigation model and ephemeris data 

- Perform band to band geometric calibration 

- Perform radiometric calibration of new sensor (relative to old one) 

General Characteristics 

• Use database of Ground Control Points (GCP) or Chips 

• Normalized Cross-Correlation (NCC) is the most common similarity measure 

• Digital Elevation Model (DEM) is rarely integrated in the registration process 

• Cloud masking usually integrated 

• Errors in the [0.15-0.5] range 
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Orthogonal Wavelet Image Registration 


Columns 




Orthogonal Wavelet Image Registration 

Rotation and Translation Invariance Issues 

• Nyquist criterion, sample signal at least twice frequency of highest 
frequency component 

- in og wavelets, signal changes within or across subbands with 
subsampling 

• Study for Shift Sensitivity (Stone et al, 1999): 

- low-pass subband relatively insensitive to translation, if features 
are twice the size of wavelet filters 

- high-pass subband more sensitive but can still be used. 



/ 
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Correlate Wavelets of Two Pulses 
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Orthogonal Wavelet Image Registration 

Rotation and Translation Invariance Issues (cont.) 


Correlation Minima vs Pulse width 
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Rotation- and Translation-Invariant Pyramids 


• Simoncelli: 

- Relax critical sampling condition of 
wavelet transforms 

- Overcomplete representation by 
4k/3 (k: number of band-pass 
filters) 

• Splines: 

- Recursive anti-aliasing prefiltering 
followed by a decimation of 2 

- Only low-pass bands 


Original 
Image or 
Output of LI 



Simoncelli 

Steerable Pyramid 



Original 
Image or 
Output of LI 



Sit 



SI 


Simoncelli with 1 Band-Pass filter (k=1) 
and 4 levels of decomposition 
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Comparative Studies Using Synthetic Data 

(Reproduced from Zavorin & Le Moigne, 2005) 



Synthetic Image Generation 









Orthogonal Wavelet Studies 

(Reproduced from Le Moigne & Zavorin, 2000) 



SHIFT ERROR CONTOURS — DAUBECHIES, LARGE ROTATIONS, T x - 0 


ROTATION (LARGE) 


Shift Errors - Daubechies - Large Rotations 


Rot = 10 T x = 20 T y = 30 



Shift Errors Function of Noise - Daubechies 
Large Rotation and Translation 
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Spline and Simoncelli Pyramids Studies 

(Reproduced from Zavorin & Le Moigne, 2005) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard deviation 
converged error 

TRU-SpIC 

TRU-SimB 

TRU-SimL 

1657/7236 ~ 22.9% 
1552/7236 = 21.5% 
3693/7236 = 51% 

0.0219008 

0.0411122 

0.0214522 

0.086284 

0.141976 

0.056263 

0.148911 

0.212933 

0.109165 

Average Error for Converged Region of Test Dataset (Warp & Noise) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard deviation 
converged error 

TRU-SpIC 

2831/9801 ~ 28.9% 

0.285565 

0.300596 

0.082751 

TRU-SimB 

725/9801 « 7.4% 

0.052036 

0.065967 

0.032609 

TRU-SimL 

2918/9801 =29.8% 

0.320529 

0.331067 

0.091969 

Average Error for Converged Region of Test Dataset (Warp & PSF) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard 
deviation 
converged error 

TRU-SpIC 

1424/7236 ~ 
19.7% 

0.216392 

0.278916 

0.168494 

TRU-SimB 

1415/7236 = 
19.6% 

0.164233 

0.252788 

0.213441 

TRU-SimL 

4038/7236 « 
55.8% 

0.243106 

0.289142 

0.142683 


Average Error for Converged Region of Test Dataset (Warp & PSF & Noise) 
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Spline and Simoncelli Pyramids Studies 

(Reproduced from Zavorin & Le Moigne, 2005) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard deviation 
converged error 

TRU-SpIC 

TRU-SimB 

TRU-SimL 

1657/7236 ~ 22.9% 
1552/7236 = 21.5% 
3693/7236 = 51% 

0.0219008 

0.0411122 

0.0214522 

0.086284 

0.141976 

0.056263 

0.148911 

0.212933 

0.109165 

Average Error for Converged Region of Test Dataset (Warp & Noise) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard deviation 
converged error 

TRU-SpIC 

2831/9801 ~ 28.9% 

0.285565 

0.300596 

0.082751 

TRU-SimB 

725/9801 « 7.4% 

0.052036 

0.065967 

0.032609 

TRU-SimL 

2918/9801 =29.8% 

0.320529 

0.331067 

0.091969 

Average Error for Converged Region of Test Dataset (Warp & PSF) 



Number of 
converged 

Median converged 
error 

Mean converged 
error 

Standard 
deviation 
converged error 

TRU-SpIC 

1424/7236 ~ 
19.7% 

0.216392 

0.278916 

0.168494 

TRU-SimB 

1415/7236 = 
19.6% 

0.164233 

0.252788 

0.213441 

TRU-SimL 

4038/7236 « 
55.8% 

0.243106 

0.289142 

0.142683 


Average Error for Converged Region of Test Dataset (Warp & PSF & Noise) 
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A Framework for the Analysis 
of Various Image Registration Components 



Strategy 







Web-based Image Registration 
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Web-based Image Registration 
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Algorithm Testing Using 
Landsat-TM Multitemporal Data 
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Algorithm Testing Using 
Multisensor Data (ETM, IKONOS and MODIS) 

Red and NIR Bands; 30m- 4m - 250 and 500 m respectively 




ETM+ 


IKONOS 
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Framework Testing Using Synthetic Datasets 

Marquart-Levenberg Optimization Using L2-Norm and Mutual Information 

(Reproduced from Zavorin & Le Moigne, 2005) 


TRiMspicaa* TRu-sm*9i% TBu-swrteas 



Contour Plot “SameRadNoisy” Dataset 
Optimization Using L2-Norm 
Threshold of 0.5 


TRUMI-SptC 51% TFlUMI-SimB 85% TRUMI-Sml 75% 



Contour Plot “SameRadNoisy” Dataset 
Optimization Using Mutual Information 
Threshold of 0.5 
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Framework Testing Using Synthetic Datasets 

Stocchastic Gradient Optimization Using L2-Norm and Mutual Information 

(Reproduced from Zavorin & Le Moigne, 2005) 


SPSA-SplC 63% $PSA-Sim6 72% SPSA-SimL 75% 



Contour Plot “SameRad Noisy” Dataset 
Stocchastic Gradient and Mutual Information 
Threshold of 0.5 



Contour Plot “SameRadNoisy” Dataset 
Fast Fourier Correlation 
Threshold of 0.5 
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Multitemporal Datasets 

Robust Feature Matching Using Simoncelli Band-Pass Features 

(Reproduced from Netanyahu et al, 2004) 
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Results of Multitemporal Registration 
Using Landsat-TM Data over DC/Baltimore Area 



RFM REGISTRATION 

MANUAL GROUND 
TRUTH 

ABSOLUTE ERROR 


Scene 

Q 

T 

1 X 

T 

V 

Q 

T 

* X 

T 

ly 

|DQ| 

IDTJ 

IDTJ 

990804 

0.009 

0.36 

3.13 

0.002 

0.04 

3.86 

0.011 

0.40 

0.73 

991108 

0.000 

LOO 

13.00 

0.002 

1.20 

13.53 

0.002 

0.20 

0.53 

000228 

0.005 

0.88 

-2.32 

0.008 

1.26 

2.44 

0.003 

0.38 

0.12 

000822 

0.002 

0.41 

9.22 

0.011 

0.35 

9.78 

0.013 

0.06 

0.56 


Results of Multitemporal Registration 
Using Landsat-TM Data over Virginia Area 
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Multisensor Datasets 

All Algorithm Comparison 

(Reproduced from Le Moigne et al, 2001) 
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Results of Multisensor Registration 
Using ETM, IKONOS and MODIS Data over Konza Agricultural Area 


Similar Tests performed on: 

- Urban Area (USDA site; Greenbelt, MD) 

- Coastal Area (VA Coast) 

- Agricultural Area (Cascades Site, CO) 

- Mountainous Area (Konza Prairie, Kansas) 

Consistency studies show between 0.125 and 0.25 pixel errors 
using circular registrations of IKONOS NIR and Red data 


Additional studies performed on EOI -Hyperion data 
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Fusion of Remotely Sensed Data 


• Data Fusion 

- Use multi-source data of different natures to increase quality of 
information contained in data (Pohl and Genderen, 1998) 

- A process dealing with association, correlation, and combination 
of data and information from single and multiple sources to 
achieve refined position and identity estimates, and complete 
and timely assessments of situations and threats, and their 
significance (Hall and Llinas, 2001). 

• Image Fusion 

- Data are images 

- General Objectives: 

» Image sharpening 

» Improving registration/classification accuracy 
» Temporal change detection 
» Feature enhancement 

- Example Application 

» Invasive Species Forecasting System 
» Objective 

» Improvement of classification accuracy 

» Tamarisk, Leafy Spurge, Cheat grass, Russian olive, etc. 

» Feature enhancement 


1 


Image Fusion Methods 


Principal Component Analysis, PCA 

- Input 

• Multivariate data set of inter-correlated variables 

- Output 

• Data set of new uncorrelated linear combinations of the 
original variable 


Wavelet-based Fusion 

- Use of Different Subbands in Reconstruction 

Cokriging 


Image Fusion Methods 

Wavelet-Based Image Fusion 


High Spatial 
Resolution 


Data 




Low Spatial 
Resolution 
Dat a — I 








Decomposition 


> 


FUSED 

DATA 


■ 
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Improved 

Resolu t ions 


Reconstruction 






Image Fusion Methods 

Cokriging 


Interpolation Method 

- Geo-statistics, mining, and petroleum engineering 
applications (pioneered by Danie Krige, 1951) 

- Generalized version of kriging (B.L.U.E): 

• Best: aims to minimize variance of the errors 

• Linear: estimates are weighted linear combination of the available data 

• Unbiased: tries to have mean residual, or error, equal to zero. 

• Estimator 

Interpolation using more that one type of variable to 
estimate an unknown value at a particular location 


Goal of cokriging is to minimize variance of error 
subject to some constraints (to ensure unbiasedness 
of our estimate) 


Image Fusion Experiments 

Using Principal Component Analysis 

(Reproduced from Memarsadeghi et al, 2005) 

Input 

- 9 bands of ALI 

- 140 bands of Hyperion (calibrated and not corrupted bands) 

- Stack of both ALI and Hyperion bands above 

Output 

- Same number of PCs as input bands 

- Select PCs containing 99% of information 



ALI PCs 1,2,3 



ALI PCS 1,2,3 
clustering 



Hyperion 
First 7 PCs 
clustering 


ALI- Hyp First 9 PCs 
clustering 


ALI V 

Hyp V 

Fused V 

143,98 

137,64 

180,10 





Image Fusion Experiments 

Using Wavelet-Based Fusion 

(Reproduced from Memarsadeghi et al, 2005) 

Fuse each multispectral band of ALI with one band of Hyperion 

• For each of 9 ALI bands 

• Select a Hyperion band within the wavelength range of corresponding ALI band which is 

» closest to the center of ALI’s wavelength range (experiment 1 ) 

» least correlated to the corresponding ALI band (experiment 2) 

• Clustering of fusion result of 9 bands of ALI with 9 bands of Hyperion 

• Fusion: 4 Levels of Decomposition, Daubechies Filter of size 2 


Hyperionl 


Fused: First test 


Fused: 2nd test 


- Experiment 1 Variances: ALI: 179.73; Hyperion: 159.96; Fused: 195.27 3 

- Experiment 2 Variances: ALI: 179.73; Hyperion: 165.34; Fused: 173.77 


Image Fusion Experiments 

Using Cokriging 

(Reproduced from Memarsadeghi et al, 2006) 



Landsat-TM 

Multispectral Bands 2, 3, 4 
(30m resolution) 



Landsat-TM Panchromatic 
(15m resolution) 


Image Fusion Experiments 

Using Cokriging (conk) 

(Reproduced from Memarsadeghi et al, 2006) 


Landsat-7 Multispectral 
Bands 2,3 and 4 
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Spectral Resolution 
1 pixel of an MS band 
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p2 

? 

x3 y3 

p3 
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x4 y4 

p4 

? 


PAN 


Landsat-7 Pan-Sharpened MS Bands 2,3 and 4 
Through Cokriging with Pan Band 8 



Results: 

• Correlation: Wavelet: 0.86; PCA: 0.91; Cokriging: 0.92 

• Entropy: Wavelet: 3.44; PCA: 3.87; Cokriging: 3.92 
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ACRONYMS 


• AVHRR: Advanced Very High Resolution Radiometer 

• CP or GCP: Control Point or Ground Control Point 

• GSHHS: Global Self-consistent Hierarchical High-resolution Shoreline 

• GOES: Geostationary Operational Environmental Satellite 

• GSFC: Goddard Space Flight Center 

• MISR: Multiangle Imaging SpectroRadiometer 

• MODIS: MODerate resolution Imaging Spectrometer 

• NDVI: Normalized Difference Vegetation Index 

• SeaWiFS: Sea-viewing Wide Field-of-view Sensor 

• SPOT: Satellite Pour I’Observation de la Terre 

• WSU: Wright State University 
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Image Fusion by Image 

Compositing 

Arthur Goshtasby 
Wright State University 
Image Registration and Fusion Systems 


Image fusion 


Combining relevant information in two or 
more images of a scene into a single highly 
informative image. 

Fusion methods: 

- Low-level (pixel-based) 

- Mid-level (region-based) 

- High-level (description-based) 


A. Goshtasby 


Image compositing 


Cutting and smoothly stitching together parts 
of different images of a scene into a single 
highly informative image. 



Image 1 Image 2 Combined image 


A. Goshtasby 



Image fusion as an image blending 

problem 



A. Goshtasby 





1 . For each local area in the image 
domain, find the image 
providing the most information. 

2. Cut out the 


Image domain 


area out of 
the most 
informative 
image. 



A. Goshtasby 
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3. Blend the local areas to create the fused 


image. 



A. Goshtasby 



Selecting the best-block image 


Highest contrast 
Highest entropy 

255 

£g=Z# 1 o S(a) 

1=0 

E c = E r +E g + E b 


Gray-scale 

Color 


A. Goshtasby 


The blending function 


Assuming (x u yj) represents the coordinates of 
the center of the zth block, use 


WXx,y) = 


GXx,y) 


as the blending function, where 

G^J ^-xS+Jy-y.y 

<7 is adjusted to the size of the local area. 


A. Goshtasby 
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Image blending algorithm 


Given images {I k (x,y)\ k= and assuming 

1 ) the image domain is subdivided into N 
blocks, and 

2) for block i, the most informative image is 

hkfay)’ 

the fused image is computed from: 

o(x, y)=Yj w i (*> yVa ( x > y) 

i = 1 


A. Goshtasby 


Characteristics of the algorithm 


• Parameter a is adjusted to block size, and 
block size is chosen based on image size 
and/or image resolution. 

• The blending process does not create 
information in output that does not exist in 
input. 

• The method is versatile, it can be used to fuse 
various types of images. 


A. Goshtasby 
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Application 1: 

Fusion of multi-exposure images 

Optimizing criterion: Entropy 


Images courtesy of Paul Debevec, USC 


A. Goshtasby 
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Optimizing criterion: Entropy 



Original images courtesy of Shree Nayar, Columbia University 


A. Goshtasby 
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Optimizing criterion: Entropy 



Original images courtesy of Max Lyons 


A. Goshtasby 
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Application 2: 

Fusion of multi-focus images 


Original 

images 



m>J!) 


5 // 


Fused 

image 



A. Goshtasby 


Optimizing criterion: 
Contrast 

Original images courtesy of Image 
Fusion Systems Research 


Optimizing criterion: Contrast 


Original 

images 


Fused 

image 




A. Goshtasby 


Over55.QOO precise definitions, 
including many new' words 

Style and usage guides for school and business 

Up-to-date* biographic and geographic entries 

Latest computer and science terms 
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Concluding remarks 


If the type of information to be maximized is known, 
this method is guaranteed to produce an image with 
the highest information content at the desired 
resolution. 

It can be adapted to image resolution by making 
block size a function of image resolution. 

By changing the optimization criterion, the method 
can be used to fuse various types of images. 


A. Goshtasby 
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Concluding Remarks 


Arthur Ardeshir Goshtasby (WSU) 

& 

Jacqueline Lemoigne (NASA) 


Summary 


Tools and methods for image registration were 
reviewed. 

Methods for the registration of remotely sensed data 
at NASA were discussed. 

Image fusion techniques were reviewed. 

Challenges in registration of remotely sensed data 
were discussed. 

Examples of image registration and image fusion 
were given. 


Conclusions 


A single method cannot register or fuse all 
image types. 

A custom-made method is require for a 
particular class of images. 

As the intensity and geometric difference 
between images increase, the registration 
process becomes more complex. 

Image fusion necessitates image registration. 


Open problems 


An image registration method is needed that can 
automatically adapt to intensity and geometric 
differences between images. 

Methods that can verify the correctness of a 
registration are needed. 

Two main steps in image registration still require 
further research. They are: 

- Reliable landmark correspondence methods. 

- 2) Image warping models that can adapt to the density and 
organization of the landmark. 


Image registration examples 


Manual registration 
Semi-automatic registration 
Automatic registration 

- Rigid registration 

- Similarity registration 

- Affine registration 

- Projective registration 

- Nonlinear registration 


Thank you ! 


